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The exchange-correlation energy in Kohn-Sham density functional theory is expressed as a func- 
tional of the electronic density and the Kohn-Sham orbitals. An alternative to Kohn-Sham theory 
is to express the energy as a functional of the reduced first-order density matrix or equivalently the 
natural orbitals. In the former approach the unknown part of the functional contains both a kinetic 
and a potential contribution whereas in the latter approach it contains only a potential energy and 
consequently has simpler scaling properties. We present an approximate, simple and parameter-free 
functional of the natural orbitals, based solely on scaling arguments and the near satisfaction of a 
sum rule. Our tests on atoms show that it yields on average more accurate energies and charge den- 
sities than the Hartree Fock method, the local density approximation and the generalized gradient 
approximations . 

FAGS numbers: 71.15.-m, 71. 15. Mb 



The solution of the quantum mechanical many-electron 
problem is one of the central problems of physics. A great 
number of schemes that approximate the intractable 
many-electron Schrodinger equation have been devised to 
attack this problem. Most of them map the many-body 
problem to a self-consistent one-particle problem. Proba- 
bly the most popular method at present is Density Func- 
tional Theory (DFT) especially when employed with 
the Generalized Gradient Approximation (GGA) 
for the exchange-correlation energy. DFT is based on 
the Hohenberg-Kohn theorem which asserts that the 
electronic charge density completely determines a many- 
electron system and that in particular the total energy is 
a functional of the charge density. Attempts to construct 
such a functional for the total energy have not been very 
successful because of the strong non-locality of the ki- 
netic energy term. The Kohn-Sham scheme where the 
main part of the kinetic energy, the single particle kinetic 
energy, is calculated by solving one-particle Schrodinger 
equations circumvented this problem. The difference of 
the onc-particle kinetic energy and the many-body ki- 
netic energy is a component of the unknown exchange- 
correlation functional. The exchange-correlation func- 
tional is thus a sum of a kinetic energy contribution and 
a potential energy contribution and partly for this rea- 
son it does not scale homogeneously |^] under a uniform 
spatial scaling of the charge density. 

It has been known for a long time, that one can also 
construct a total energy functional using the first-order 
reduced density matrix. Several discussions of the ex- 
istence and the properties of such a functional can be 
found in the literature [|7|j9|-[Tl[| . However in spite of the 
enthusiasm expressed towards this approach in the early 
papers, no explicit functional has ever been constructed 
and tested on real physical systems. An important ad- 



vantage of this approach is that one employs an exact 
expression for the many-body kinetic energy. Only the 
small non Hartree-Fock-like part of the electronic repul- 
sion is an unknown functional We propose in this 
paper an explicit form of such a functional in terms of 
the natural orbitals. The high accuracy of this Natural 
Orbital Functional Theory (NOFT) is then established 
by applying it to several atoms and ions. 

Let us first briefly review some basic facts about re- 
duced density matrices . If ^' is an arbitrary trial 
wave function of an A^-electron system, the first and sec- 
ond order reduced density matrices, 71 and 72 are 



7i(x;,xi) = TV J ... J *(x;,X2,...,XAr) 



*(xi,X2, ...,XAr) (ix2...dxAr, 

iV(iV- 1) 



N(N - 1) /■ f 
72(x'i,x^;xi,X2) = / ... / *(x;,x^,X3,. 



^'(xi,X2, ,X3, ...,XAr) dX3...dx 



(1) 

,XAr) 
(2) 



The variables contain both the position coordinates 
Vi, as well as the spin coordinate s^. The integration sign 
stands for a combined integration of the spatial coordi- 
nates and summation of the discrete spin part. 

The electronic charge density p{r) is obtained from the 
diagonal part of the first-order reduced density matrix. 



p(xi) = 7i(xi,xi); p(ri) = ^p(xi). 



(3) 



The natural orbitals (pi are the eigenfunctions of the 
first-order reduced density matrix with eigenvalues rii. 



y"7i(x'i,xi)0i(xi)dxi =ni<j)i{x[) 



(4) 



The natural spin-orbitals and occupation numbers Ui 
specify the reduced first-order density matrix completely. 
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The total energy can be written in terms of the natu- 
ral orbitals and the diagonal elements of the second order 
reduced density matrix, 



E : 



Cr(xi,X2) = 72(xi,X2;Xi,X2), 

-^Yl^i J 0j(x)V2(?!),(x)dx 



(5) 



(6) 



ri - r2 



(ixi(ix2 



In order to construct a natural orbital functional, it 
remains to find an approximation for a in terms of the 
natural orbitals and occupation numbers. In the follow- 
ing, we assume the standard case of a Hamiltonian that 
is not spin dependent. Each natural orbital can then be 
chosen to be either purely spin up or spin down and can 
be labeled by an orbital index i and a spin index s^. 

The approximate a we propose has the following form: 



^[M,W]=E:,^0?(ri)0,^(r2) 



(7) 



-Eij ^^^ '?^.,^,0»(ri)0j(ri)0,(r2)(/)j(r2) . 

The primes indicate that the i — j terms are omitted. To 
find the ground state, we minimize the functional with 
respect to both the natural orbitals and the occupation 
numbers, under the constraint that the natural orbitals 
be orthogonal The functional derivatives are 



-|^ = -^V20.(r)+n,y(r)<^,(r) 



(8) 



r — r' 



|r — r'l 



^,=-\j Uv)V^U^)dv + j V{v) 4>l{v) dv (9) 
</.?(r')^?(r) 



drdr' 



</>,(r')</>,(r')0.(r)<^,(r) 



drdr' 



In principle an infinite number of natural orbitals must 
be included. For the systems studied in Table | at 
most 38 orbitals were needed to obtain good convergence. 
The occupation numbers of the core natural orbitals are 
restricted to be unity, while the remaining occupation 
numbers are allowed to vary freely and are found to 



lie always between zero and one, which is a necessary 
and sufficient condition for the density matrix to be N- 
representable [ p^ . 

We now discuss the properties of this functional. 

Homogeneous scaling of exchange-correlation energy: 
The exact exchange-correlation energy in first-order den- 
sity matrix functional theory differs from the exact 
exchange-correlation energy in density functional theory 
amd scales homogeneously under a uniform scaling 
of the density matrix. The exchange-correlation energy, 
deduced from Eqs. (^) and (0), exhibits this property. 

No orbital self-interactions: 
In the case where one has fractional occupation numbers 
one has to distinguish between orbital self-interactions 
and electron self-interactions. Our functional is free of or- 
bital self-interactions because the sum in Eq. |^ excludes 
terms with i — j, but it is not perfectly electron self- 
interaction free. The total energy for H is therefore not 
correct (Table |). The functional has however a much 
better cancellation of electron self-interactions than den- 
sity functionals, as can be seen from the fact that nega- 
tive ions are stable (Table |) . In contrast LDA and GGA 
bind only a fraction of an additional electron. 

Sum rule for second order reduced density: 
The density and the number of electron pairs are ob- 
tained by integrating the exact second order reduced den- 
sity matrix. 



/ cr(ri,r2) dr2 = p{ri 



cr(ri,r2) dridr2 



N{N- 1) 



(10) 
(11) 



Our approximation for the second order reduced density 
matrix would satisfy these equations if the sums in Eq. |^ 
also included the i = j terms. We omit these terms be- 
cause we find that an exact cancellation of the orbital 
self-interactions is more important than an exact fulfill- 
ment of the sum rules in Eqs. ( p^ ) and (pi]). The sum 
rules are violated only by terms of the order of ni(l — n^), 
which for most systems are small since all the occupation 
numbers are close to either zero or one. 

Hartree Fock as limiting case: 
The functional coincides with the Hartree Fock (HE) 
functional if one imposes the additional constraint, that 
the occupations numbers all be 1 or 0. 

No dissociation problems: 
Even though the functional contains terms which are 
similar to the HE functional, it should not suffer from 
some well established deficiencies of the spin restricted 
HE functional such as the dissociation problem of the H2 
molecule. As one separates the two H atoms, the large 
occupation numbers in the up- and down-spin ag molec- 
ular orbital get redistributed to the up-spin Is atomic 
orbital on one atom and the down-spin Is atomic orbital 
on the other. In the infinitely separated limit each atom 
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has non-zero occupation numbers in either only the up- 
spin or only the down-spin orbitals. Consequently the 
energy is the sum of the energies of the individual atoms. 

Transition states: 
In molecular calculations the effect of this functional 
is expected to be particularly significant for transitions 
states, which are poorly described by LDA and HF. At 
transition states more than one determinant is needed for 
an adequate description, and releasing the HF constraint 
of integer occupation numbers is therefore important. 

Orbital-dependent "potentials" : 
The weakly-occupied natural orbitals are localized in the 
same region of space as the highest strongly-occupied 
natural orbitals. This is in contrast to the unoccu- 
pied Kohn-Sham and Hartree-Fock orbitals which have 
a larger extent than the occupied ones. The manner in 
which this comes about can be seen from Eq. ^ which 
has an orbital-dependent "potential" . One term in the 
potential goes as ^Jnl - an enhancement by a factor of 
l/^/rTi relative to Hartree-Fock - which has the conse- 
quence that weakly-occupied natural orbitals see a more 
strongly negative potential than do the strongly-occupied 
orbitals, thereby helping to localize the weakly-occupied 
natural orbitals. 

Chemical potential: 
All natural orbitals with fractional occupation rii share 
the same chemical potential , /i = ^ . 

Discontinuity of the exchange-correlation potential: 
As one adds fractions of an electron, one finds, at occu- 
pation numbers close to integers, a rapid change in the 
effective potential felt by all the electrons, which is due to 
the jump in the chemical potential. This quasi discontin- 
uous effect might mimic the discontinuity in the 
DFT exchange correlation potential, an effect missing in 
the LDA and GGA functionals. 

Correct description of correlations of different origin: 
In a 1/Z expansion of the energy, the correlation energy 
of the two-electron series can be described by nondegen- 
erate perturbation theory while the four-electron series 
requires degenerate perturbation theory. Consequently 
the correlation energy of the two-electron series tends to 
a constant with increasing Z, whereas it increases lin- 
early in the four-electron case. Both trends are correctly 
captured by the NOFT functional as shown in Table ||. 
Any GGA functional can at best describe only one of the 
trends. 

Correct qualitative behavior of natural occupation 
numbers: 



As seen from Table HI, the NOFT occupation numbers 



may differ considerably from the ones obtained from con- 
figuration interaction calculations, but the main trends 
are correctly reproduced. In particular, the trend in the 
occupation numbers of the strongly occupied Is orbitals, 
going from He to H~ is correct. 

Accurate results: 
In Table |, we give a compilation of the errors in the 



total energy Ai? and the errors in the charge densities 
Ap. The charge density errors are defined by (Ap)^ = 
/(/Ocx(r) — /o(r))^(ir, with the "exact" charge densities 
Pcx obtained from accurate quantum Monte Carlo calcu- 
lations Both total energies and charge densities are 
improved on average compared to HF and DFT calcula- 
tions. In particular the improvements over the HF densi- 
ties are impressive since they are known to be rather ac- 
curate. The GGA schemes yield improved total energies 
compared to both LDA and HF while the GGA densities 
are better than those from LDA but not as good as those 
from HF. In the case of C, the error in the spherically av- 
eraged charge density is quoted. The "exact" total ener- 
gies were obtained from Ref . . The LDA energies and 
densities were obtained by a standard spherical atomic 
program. As a representative of a GGA functional we 
have chosen the recent PBE functional. All the HF 
and NOFT calculations were done with a non-spherical 
atomic program developed by the authors. All calcula- 
tions were done in a spin restricted scheme. In the case of 
C the correct non-spherical '^P ground state was chosen. 
The QCISD configuration interaction calculations were 
done with the Gaussian 94 software package using 
an accurate 6-311G-f -|-G(3df,2p) basis set. Since we do 
no molecular calculations, we monitor a third quantity 
the transferability error At, to make predictions about 
the behavior of this scheme in molecular and solid state 
calculations. Molecular geometries are determined via 
the Hellmann-Feynman theorem by the charge densities 
in the valence region. The external potential in the va- 
lence region is modified in a molecule compared to the 
atomic case. We simulated this modification by adding 
a confining parabolic potential to the atom. The change 
in the total energy due to the variation of this parabolic 
potential is again given by the Hellmann-Feynman theo- 
rem and we define the transferability error Ar therefore 
as At = / (p(r) - Pox(r)) r'^dv. 

In conclusion, we have made a first attempt at con- 
structing an approximate total energy functional of the 
first-order reduced density matrix. We have listed and 
discussed the properties that make it superior to the HF 
and approximate DFT functionals and have also shown 
that it yields better energies and densities than HF and 
current DFT schemes. The high accuracy of quantities 
related to the charge density leads one to expect that 
this new functional will give accurate molecular geome- 
tries as well as accurate energy differences between dif- 
ferent geometric configurations. In view of the fact that 
the functional is parameter free and based on a few sim- 
ple considerations, we think this to be a remarkable suc- 
cess. It is likely that it will be possible to construct even 
better functionals along these lines. The essential point 
in this work is that we have used natural orbitals in- 
stead of Kohn Sham orbitals. We believe that this is es- 
sential to obtain accurate densities and kinetic energies. 
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With the exception of the complication of an orbital- 
dependent potential, the computational procedure neces- 
sary to solve the NOFT equations is analogous to other 
self-consistent one-particle schemes and thus computa- 
tionally much cheaper than quantum chemistry methods 
based on configuration interaction related schemes. 

We thank M. Levy, O. Gunnarson, J. Hutter, M. 
Teter, K. Maschke, A. Savin and B. Fahrid for interest- 
ing discussions and for suggesting references. Mike Harris 
kindly provided a subroutine for angular grid generation. 
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TABLE I. Comparison of the errors of the quantities de- 
scribed in the text. Energies are in Hartree atomic units. No 
data are available (NA) for the non-spherical PBE ground 
state of C. The large errors in Ap and the infinite errors in 
At for the H~ ion in LDA and PBE come from the fact that 
they bind only a fraction of the additional electron. 





H 




He 


Li 


Be 


C 


Ne 




Energy 


- E 


.5 


.5278 


2.9037 


7.4781 


14.6674 


37.8450 


128.9376 




LDA 


AE 


2.e-2 


6.e-3 


7.e-2 


l.e-1 


2.e-l 


4.e-l 


7.e-l 




l.e-3 


6.e0 


8.e-3 


2.e-2 


2.e-2 


5.e-2 


2.e-l 


At 


4.e-l 


oo 


2.e-l 


-7.e-l 


2.e-2 


4.e-l 


3.e-l 




PBE 


AE 


8.e-5 


2.e-3 


l.e-2 


2.e-2 


4.e-2 


NA 


7.e-2 


(Ap)^ 


2.e-4 


6.e0 


l.e-3 


3.e-3 


3.e-3 


NA 


l.e-2 


At 


2.e-l 


oo 


l.e-1 


-l.eO 


5.e-l 


NA 


3.e-l 




HE 


AE 


0. 


4.e-2 


4.e-2 


5.e-2 


9.e-2 


2.e-l 


4.e-l 


(Ap)' 


0. 


l.e-3 


l.e-4 


7.e-5 


8.e-4 


5.e-4 


6.e-3 


At 


0. 


-5.e0 


-2.e-2 


3.e-l 


l.eO 


-6.e-2 


-2.e-l 




NOFT 


AE 


-2.e-2 


l.e-2 


6.e-3 


-l.e-3 


-2.e-2 


3.e-2 


5.e-2 


(Ap)' 


3.e-5 


4.e-4 


l.e-5 


2.e-4 


6.e-4 


7.e-4 


4.e-4 


At 


-2.e-2 


l.e2 


-l.e-2 


-5.e-l 


6.e-l 


5.e-2 


-5.e-2 




QCISD 


AE 


2.e-3 


8.e-3 


8.e-3 


5.e-2 


5.e-2 


7.e-2 


l.e-1 



TABLE II. Correlation energies, in Hartrees, for 
the 2- and 4-electron series. The exact values of 
-E^'^ = E^^ - taken from Ref @, are compared 

to SHF_^NOPT^ 



2 electron 4 electron 



z 






Z 






1 


.040 


.031 


4 


.094 


.110 


2 


.042 


.036 


6 


.126 


.141 


4 


.044 


.040 


8 


.154 


.171 


6 


.045 


.042 


10 


.180 


.200 



TABLE III. Occupation numbers for the 2-electron se- 
ries. Columns labeled 'E' are the almost exact numbers of 
Kutzelnigg [M. Entries smaller than le-5 were set to zero. 



nl 


E:Z=1 


Z=l 


E:Z=2 


Z=2 


E:Z=4 


Z=4 


E:Z=6 


Z=6 


Is 


.9646 


.9666 


.9930 


.9943 


.9984 


.9985 


.9993 


.9993 


2s 


.24e-l 


.lOe-1 


.32e-2 


.22e-2 


.63e-3 


.53e-3 


.25e-3 


.24e-3 


2p 


.lle-1 


.28e-2 


.34e-2 


.92e-3 


.89e-3 


.27e-3 


.39e-3 


.13e-3 


3s 


.73e-4 


.lOe-1 


.36e-4 


.lle-3 





.19e-4 








3p 


.15e-3 


.47e-3 


.79e-4 


.77e-4 


.23e-4 


.16e-4 








4s 





.27e-3 




















3d 


.37e-3 


.30e-3 


.15e-3 


.61e-4 


.47e-4 


.13e-4 


.21e-4 





4p 


.12e-4 


.lle-3 




















5s 





.81e-4 




















4d 


.33e-4 


.57e-4 


.14e-4 

















5p 





.91e-4 




















6s 





.12e-4 




















6p 



























